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Abstract. The Nose-Hoover thermostat is a deterministic dynamical system designed for comput- 
ing phase space integrals for the canonical Gibbs distribution. Newton's equations are modified by 
coupling an additional reservoir variable to the physical variables. The correct sampling of the phase 
space according to the Gibbs measure is dependent on the Nose-Hoover dynamics being ergodic. 
Hoover presented numerical experiments that show the Nose-Hoover dynamics to be non-ergodic 
when applied to the harmonic oscillator. In this article, we prove that the Nose-Hoover thermostat 
does not give an ergodic dynamics for the one-dimensional harmonic oscillator when the "mass" of 
the reservoir is large. Our proof of non-ergodicity uses KAM theory to demonstrate the existence 
of invariant tori for the Nose-Hoover dynamical system that separate phase space into invariant 
regions. 

We present numerical experiments motivated by our analysis that seem to show that the dynam- 
ics is not ergodic even for a moderate thermostat mass. We also give numerical experiments of the 
Nose-Hoover chain with two thermostats applied to the one- dimensional harmonic oscillator. These 
experiments seem to support the non-ergodicity of the dynamics if the masses of the reservoirs are 
large enough and are consistent with ergodicity for more moderate masses. 



1. Introduction 



Equilibrium statistical properties of molecular systems [3, 10] are given by phase space integrals 
of the form 



{A) = / A(q,p)d»(q,p), 



(1.1) 



where q = (q%, . . . , qu) £ W 1 and p = (pi, . . . ,Pm) £ K denote a set of positions qi G K n and 
momenta pi £ W 1 of M particles (n denotes the space dimension), and A(q,p) is an observable, a 
function defined over the phase space and related to the macroscopic quantity under study. If the 
molecular system is observed at fixed temperature 9, then the measure d[x is the Gibbs measure for 
the canonical ensemble [3, 10] 



dfj,(q,p) 



exp(-(3H(q,p)) 
exp (— (3H(q,p)) dqdp 



dq dp, 



(1.2) 
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where H(q,p) is the Hamiltonian of the system and is often simply of the form 

M 2 
i=l * 

where pf = Pi-pi and V(q) is the potential energy. The parameter (3 that appears in (|1.2|) is related 
to the temperature 9 by f3 = l/(ksO), where is the Boltzmann constant. In applications of 
interest, the number of particles is often very large (M > 100, 000), hence computing integrals such 
as is a challenging problem. 

Molecular dynamics can be used to compute integrals such as (jl.lj) . The method amounts to 
finding a dynamics on (q,p) which is ergodic with respect to the measure d[i. As a consequence, 
the phase-space average can be replaced by a time average 



/ A(q,p)exp(-f3H(q,p))dqdp 1 . T 

I A(q,p)dn{q,p) = - f = lim - / 

I exp(-f3H(q,p))dqdp + °° 1 Jo 



A(q(t),p(t))dt (1.3) 



over a trajectory (q(t) , p(t))t>o ■ The time average can be approximated by a formula such as 



r 



If 1 

lim -/ A(g(t),p(t))cft« lm — JjA(g*p/), 



where (qi,pe)e>i is a numerical solution of the chosen dynamics. 

To compute phase space integrals in the canonical ensemble, several deterministic dynamics have 
been proposed, such as the Nose [12], the Nose-Hoover [4], and the Nose- Hoover chain dynamics [8]. 
More recently, the Nose-Poincare dynamics [1] and the Reversible Multiple Thermostat method [7] 
have been proposed. Stochastic dynamics (such as the Langevin equation) can also be considered, 
although we will not discuss them in the following. 

The ergodicity condition has not been rigorously proven for any of the deterministic methods 
mentioned above. In fact, there is numerical evidence that shows that the Nose and the Nose- 
Hoover methods are not ergodic for some systems [4,8,14], including the one-dimensional harmonic 
oscillator. This article further explores non-ergodic behavior of the Nose-Hoover dynamics in this 
simple example. 

In Section |5J we present the Nose-Hoover equations and we recall some of their properties. In 
Section |HJ we prove that the Nose-Hoover thermostat does not give an ergodic dynamics for the 
Gibbs measure (|1.2[) for the one-dimensional harmonic oscillator when the "mass" of the reservoir 
is large. Our method is to apply KAM theory, and more specifically Moser's invariant curve 
theorem [13], to the Poincare return map and to thus demonstrate the existence of invariant tori 
that separate phase space into invariant regions. Finally, in the last section, we present some 
numerical experiments with a Nose-Hoover chain of two thermostats. For large reservoir masses, 
results show that the dynamics seems to be non ergodic. For moderate reservoir masses, they are 
consistent with ergodicity. 
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2. The Nose-Hoover Thermostat 



The Nose-Hoover dynamical system [4] is given by 

dqj _ Pi 
dt rrii ' 




1 2 



(2.1) 



where the phase space is described by the physical positions q = (q±, . . . , qu ) £ M. nM and momenta 
p = (pi, . . . ,pm) £ W 1 and an additional variable, £, which can be considered as the momentum 
of the thermostat. The constant Q, which is a parameter of the method, represents the mass of the 
reservoir and describes the strength of the coupling of the reservoir to the physical system. Let us 
note that, usually, a second additional variable is introduced [4]. This variable can be considered 
as the position of the thermostat. Since it is decoupled from all the other variables, we ignore it in 
the following. 

We recall that for the canonical Gibbs measure d[i given by (|1.2|) we have 



for almost all initial conditions for any dynamics on (q,p) which is ergodic with respect to the 
measure d\x. Thus, the right hand side of the Nose-Hoover dynamical equation for % is equal 
to twice the difference between the instantaneous kinetic energy of the physical system and the 
time-averaged kinetic energy of the physical system at temperature = kg 1 ft -1 with nM degrees 
of freedom. Hence, we see that if the kinetic energy of the physical system is too high for a 
sufficiently long time, then the "thermostat" added to the physical momentum equations applies 
a frictional force to damp the system. If the kinetic energy of the physical system is too low for a 
sufficiently long time, then the "thermostat" added to the physical momentum equations applies 
an "anti-frictional" force to add kinetic energy to the system. 

The Nose-Hoover system is not a Hamiltonian system and the Lebesgue measure dq dp dt; is not 
an invariant measure for the dynamics. Instead, the equations preserve a different measure which 
we will now describe. We recall that invariant measures p(z)dz for a general dynamical system 

are determined by the equilibrium equation 




so we have that 




For the Nose-Hoover system Q2.1J1 with z 




(2.2) 
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satisfies the condition (|2.2|) . so a normalized invariant measure for the Nose- Hoover system (|2.1(l is 
given by 

exp (-(3 H(q,p) + ^ 



/exp (-0 



dq dp d^ 



dq dpd^. 



(2.3) 



We recall that the flow of equations (|2.1j) is ergodic (or metrically indecomposable) with respect to 
the measure d/^NH if the phase space, M 3 , cannot be decomposed into two complementary invariant 
subsets, each with positive measure [5]. If the flow of (|2.1|) is ergodic with respect to the measure 
d^NH) then given any integrable function A G L^g^nh) the Birkhoff ergodic theorem shows that 
for almost all initial conditions 



A(q(t),p(t),£(t))dt= / A(q,p,OdfJ, N u 



lim — 

T^+oo T 

where (q(t),p(t),^(t)) is a solution of the Nose-Hoover equations (|2.1|) . 

If A = A(q,p) is an observable which depends only on the physical variables (q,p), then 

J A(q,p) cfyiNH = J A(q,p) dp, 

where dp is the Gibbs measure (|1.2|) . To see this, we observe from Fubini's Theorem that 

J A(q,p)dp N H 



(2.4) 



(2.5) 



j A(q,p)exp (-ft 




^ dqdpd^ 


J 


f exp 




Start + § 


^ dqdpd£ 



J A(q,p)exp(-f]H(q,p))dqdp ■ j 


^ exp 


(-'[ 


e- 

2Q 




J 


f exp(-(3H(q,p))dqdp- j 


f exp ^— /3 


\^ 
2Q 





(2.6) 



A(q,p) dp,. 



We thus have that if the flow of 1)2- lj) is ergodic with respect to the measure dp^n an d if the 
observable A = A(q,p) does not depend on £, then for almost all initial conditions we have from 
(TO!) that 



^im i / A(q(t),p(t))dt = / A(q,p)dp 

-^+00 i J J 



J A(q,p)exp(-(3H(q,p))dqdp 



exp(— (3H(q,p)) dq dp 



(2.7) 



where (q(t) , p(t) , £(t)) is a solution of the Nose-Hoover equations (|2.1|) . 

This derivation has been made under the assumption that the Nose-Hoover dynamics is ergodic 
with respect to the measure c£//nh- Numerical experiments show that the equality 1)2. 7j) does not 
always hold, even for long times T (in the limit of computationally reachable times). In particular, 
it is observed in [4, 8, 14] that if the system under consideration is a one-dimensional harmonic 
oscillator, that is, if n = M = 1, mi = 1, and V(q) = \q 2 , then for lots of initial conditions, there 
exist c and C with < c < C such that the corresponding solution of ()2.1j) satisfies 



c< q 2 (t) +p 2 (t) < C for all t. 



(2i 
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This fact is observed for a wide range of values of Q, including Q = 1. This behavior contradicts 
(|2.7|) . which gives in this case 

T I A(q,p)exp (-(3 ( q +P X\ dqdp 

Hm^- J A( ff (t),p(t)) eft = ^— /9 2 +P 2 \\ ' (2 ' 9) 

1 exp I — /3 I I I dqdp 



For example, if A(q,p) is a positive function whose support lies in the disk q 2 +p 2 < c, the left-hand 
side of Q2.9J1 will be zero while the right-hand side will be positive. 

These numerical experiments give an indication that the Nose-Hoover thermostatted harmonic 
oscillator is not ergodic. In the next section we will give a rigorous proof of non-ergodicity if the 
reservoir mass Q is sufficiently large. We will apply KAM theory [13] to demonstrate the existence 
of invariant tori that separate the phase space into invariant regions of positive measure. The pro- 
jections of these invariant regions to the (g,p)-plane satisfy inequalities of the form (|2.8j) and this 
explains the numerical observations. Motivated by our analysis, we show that we can find initial con- 
ditions such that the trajectory does not even sample the whole ring { (q,p) 6R 2 : c < q 2 + p 2 < C} 
for some < c < C, but only a part of it (see Figure 131 below) . 

3. Invariant Tori for the Nose-Hoover Harmonic Oscillator Dynamics 

We now write the Nose- Hoover equations in the case of a one-dimensional harmonic oscillator. 
To simplify the notation, let us assume that the particle mass is mi = 1 and the target temperature 
is such that = 1/(/cb#) = 1- In view of (|2.1|) . the system of differential equations is given by 

P = -q-e 2 &, (3.1) 

e = P 2 -i, 

where e = 1/ 1 \J~Q. 

We can introduce action-angle variables for the oscillator 

g = V / 2rcos# and p=—V2rsm6 (3.2) 

and rescale via a = e£ to get: 

6 = 1 — ea sin 6 cos 9, 

f = -2eTasm 2 6, (3.3) 
d = e(2rsin 2 6> - 1). 
These equations preserve the volume element 

dft = A(t, a) d6 dr da, 

where 

A(r,a) =eT T -° ? l 2 . 

We now make a change of variables as in the averaging method [13]. Setting 



r = f + era. sin 9 cos 9 and a = a — ef sin 9 cos 9 
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gives a new ODE of the form: 

= 1 -edsin0cos0 + O(e 2 ), 

f = - e fa + 0(e 2 ), (3.4) 
a = e (f -l)+0(e 2 ). 

The displayed terms in f , a are the averages with respect to of the corresponding terms in 
()3.3() . These equations preserve the volume element obtained by transforming dVl: 



where 
Indeed, 



A e (0,f,d) 



-r-a 2 /2 



dCl e = A e (0, f , a) d8 df da, 
(d,f,a) = e^- &2 / 2 + 0(e). 
<9(t, a) 



d(f, a) 



= e - T - a2 / 2 fl + -asin20 + — f sin 2 20 
V 2 4 

= e - f - A2 / 2 e-^ 2 sin2 9 ™ 2 8 + id sin 20 + sin 2 20^ 

= e ---« 2 /2 + ( e )_ 

In what follows, the * will be suppressed and the new variables will again be called (0, r, a). 

We will apply the KAM theory [13] to the Poincare return map P e of the plane £ = {(0,r, a) : 
= mod 2-7t} for the ODE (|3.4|) . This map preserves the area-element: 

doj € = A°(t, a) dr da, (3-5) 

with A°(r,a) = A e (0,r 5 a) = e^"" 2 / 2 . 

This can be shown as follows. Consider a small rectangle DcE with dimensions 5r5a centered 
at some point (0, TQ,ao) £ S. Form a three-dimensional tube T by following the solutions of the 
ODE (|3.4|) until they reach the plane = 2ir. One end of the tube will be the rectangle D and the 
other will be its image P e (D) under the Poincare map. 

Following T forward under the flow of the ODE for a time 5t produces a new tube T' and the 
volumes of T and T with respect to the volume element dVt e are equal. Now T and T' differ by two 
small solid "cylinders" with bases D and P e (D), respectively. Let 5 = max(|<5T|, \5a\, \St\). Then 
the volume of the cylinder over D is 

A e (0, r , a ) 5t 5a 56 + 0{5 A ) = A e (0, r , a ) 0(0, r , a ) 5t 5a 5t + 0(5 4 ), 

where is the first component of the vectorfield 1)3.4(1 . On the other hand, the volume of the 
cylinder over P e (D) is 

A e (2vr,r 1 ,a 1 ) \DP 6 (t , a )\ 5t 5a 56+0{5 A ) = A e (27r, n, c*i) 0(2tt, n,^) |DP £ (r 0) a )| <5r 5a <5t+0(<5 4 ) 

where P e (ro,ao) = (ri,ai) and where |Z?P e (ro, ao)| is the Jacobian determinant of the Poincare 
map. Note that 0(0, r, a) = 9(2ir,T,a) = 1 since this holds for the ODE 1)3. 3() and is preserved by 
the coordinate change leading to 1)3.4(1 . It then follows from letting 5 — * that 

A e (0, T ,a ) = \DP e (T ,a )\ A e (27r, n, ai) 

which proves that the Poincare map P e preserves the area element A°(r, a) dr da as claimed. 
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We will use the version of Moser's invariant curve theorem from "Lectures on Celestial Mechanics" 
by Siegel and Moser [13, sections 32-34]. That theorem starts with a real-analytic map P of the 
form 

x x = x + jy + f(x,y), 

yi = y + g{x,y), 

where / and g are periodic in x with period 2ir. In the application they have in mind, x and y are 
respectively the angle and radius in a polar coordinate system near a fixed point. If / = g = 0, then 
the map reduces to a standard twist map where the radial variable is preserved while the angular 
variable is rotated by an amount which depends on the radius. They assume a twist condition, 

7 7^ 0, so the amount of rotation really does change. The theorem shows that if / and g are small, 
then some of these invariant circles will persist. 

In order to prove this, they also assume that P satisfies the curve-intersection property. This 
means that for any simple closed curve, C, of the form 

y = ip(x) where ip(x + 2ir) = ip(x), 

we have that C n P(C) ^ 0. 

Such a curve represents a simple closed curve around the fixed point for the map before changing 
to polar coordinates. If this map preserves an area element, then such a loop cannot map completely 
inside or completely outside itself and so the curve intersection property will hold. 

To state the theorem precisely, fix an annulus a < y < b with b — a = 1 (this condition can 
always be achieved by rescaling y). Since / and g are real analytic, they extend to complex 
analytic functions on some complex neighborhood D of M. x [a, b] in C 2 . We also specify a so-called 
Diophantine condition on the rotation numbers of the unperturbed circles. Recall that if f(x) 
satisfies f(x + 2ir) = f(x) + 2ir and defines a homeomorphism of the circle M mod 2tt, then the 
rotation number u = lim n ^ 00 2^n(f n {x) — x) exists and is independent of the initial condition, x. 
For the unperturbed circles of constant y we have f n (x) — x = njy, so uo = 3^. Fix any constants 
Co > and /i > 2 and consider rotation numbers satisfying: 

>^ for all k, I e Z, I > 0. (3.7) 

It is shown in Siegel-Moser that for cq sufficiently small, the set of y which satisfy this condition 
forms a positive measure Cantor set in the interval [a, b] whose measure tends to b — a = 1 as 
Co — ► 0. Moser's theorem states that if / and g are sufficiently small, then for each such y there is 
a nearby invariant curve of P with the same rotation number. More precisely, 

Theorem 3.1. Let P be a real analytic map of the form, \3. 6\) which extends to a complex domain 
D. Assume that P is periodic in x of period 2ir and has the curve intersection property. Fix a 
Diophantine condition \3. 7[ ) and an annulus a < y < b where b — a = 1. 
Given any e > there is a 5 > such that if 

|/| + M<^7, (3-8) 

then for every yo £ [a, b] which satisfies \°J. 7| ) ; P has an invariant curve of the form y = ip{x) where 
ip(x + 2ir) = ip(x) with \ijj{x) — y$\ < e for all x. Moreover, the restriction of P to this curve has 
rotation number uj = Here \ ■ \ denotes the sup norm in the complex domain D. The constant 

8 depends on co,fi,e,D but not on 7. 

To apply this result to P e , some further coordinate changes will be needed. First, note that 

P e (r,a) =Q t (r,a)+0(e 2 ), 



\Iuj — k\ 



2vr 



NON-ERGODICITY OF THE NOSE-HOOVER THERMOSTATTED HARMONIC OSCILLATOR 



8 



where Q e is the time-2-7r advance map of the differential equation: 

f = —era, 
a = e(r — 1), 

or equivalently, the time-2-7re advance map of 

t' = —ra, 

(3-9) 

a = t — 1. 

This follows because the return time to the section S is 

T(T,a) = 2-K + 0{e). 

The ODE has an integral 

G(r,a) =t -\nr + \a 2 - I. (3.10) 

Furthermore there is an equilibrium point at (ro,«o) = (1>0) on the level set G = 0. All of the 
other level sets in the half-plane r > are simple closed curves around the equilibrium point 
(see Figure ^) . These are all invariant curves for the map Q e , and the goal is to show that the 
actual Poincare maps P e with e sufficiently small also have such invariant curves. We note that the 
existence of invariant curves for P e corresponds to the existence of invariant tori for the 3D flow of 
(|3.4|) . These tori separate the phase space into invariant 3D regions, showing that the flow is not 
ergodic. Note that for an invariant torus close to a level set G = go > we have 

r-lnr < 1 + 5-0 + 0(e)- 

Since r = + q 2 ), it follows easily that the projection of the torus to the (g,p)-plane satisfies 
bounds of the form (|2.8|) . 

It is worth noting that even without using the KAM theory, the existence of the integral G for 
the averaged system (|3.9j) shows that the convergence of ergodic averages in ()2.9j) would be slow 
for e small. Indeed, the standard averaging theory shows that the trajectories of the actual system 
remain e close to those of the averaged system on a time scale of order 1/e. The KAM theory shows 
that the trajectories which lie on the invariant tori are close to the averaged ones for all time. 
Moreover, other trajectories are trapped between the tori and so their G coordinates are prevented 
from wandering very far. 

To get the Poincare map into the form (|3.6|) . we introduce action-angle variables around the 
equilibrium point. The integral G provides a natural radial coordinate or action variable. To 
construct the corresponding angle variable, let T(g) denote the period of the periodic solutions of 
()3.9|) which correspond to the level curve G = g. We define an angular coordinate 4> to be the time 
along this orbit multiplied by 2n/T(g), taking the initial point (f> = along the r-axis to the right 
of the equilibrium point. By definition, the ODE l)3.9j) will become: 

0' = 27r/r(G), 
G' = 0, 

in these coordinates and so the time-27re advance map Q e (</>, G) = (<f>i, G%) is 

<h. = <f> + 2™/T(G), 

Gx = G. [6Ai) 
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Figure 1. Level curves of G (see (j3.10j) ), the integral for the averaged equations ()3.9|) . 



The Poincare maps P e can be viewed as 0(e 2 ) perturbations of (|3.11j) . To see this, let us first 
show that P e has a fixed point (r t ,a t ) = (1,0) + 0(e). Since Q e is the time-27re advance map of 
(|3.9jh we compute that 

Q e (r, a) = (r - 2vrera, a + 2vre(r - 1)) + 0(e 2 ), 

hence 

P e (r, a) - (r, a) = 2vrei?(e, r, a), (3.12) 

with 

i?(e, r, a) = (—to + eii(e, r, a), r — 1 + (e, r, a)) 

for some smooth functions u and v. We then apply the implicit function theorem to R(e, r, a) to 
continue the solution in e from i?(0, 1,0) = 0. We can compute 

OR , 

-(0,r,a) = 



<9(r, a) 



-a —t 
1 



so the matrix 



dR 



-(0, 1,0) is invertible. As a consequence, the implicit equation R(e,T e ,a e ) = 

o(t, a) 

defines a function e i— > (T e ,a e ) from a neighborhood of to a neighborhood of (1,0). In view of 
P-12|) . we see that (r e , a t ) is a fixed point of P t . From the equation i?(e, r e , a e ) = 0, we obtain that 
(r e ,a t ) = (1,0) + 0(e). 

After a translation of the coordinates, one can assume that all maps P e for e sufficiently small 
fix the point (1,0). Now use the same coordinates (0, G) as for the unperturbed map Q t . G is not 
an integral but one has that P t (4>-, G) = (4>i, G\) with 

1 = + 27re/T(G) + O(e 2 ), 

Gi = G + 0(e 2 ). 

The fact that P e preserves the area element du> t (see (|3.5j) ) implies that ()3.13|) satisfies the curve- 
intersection property. We also need a twist condition on the period T(G). It will be shown below 
that T'(G) > 0, that is, the period of the periodic solutions increases as we move out from the 
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equilibrium point. This implies that one can replace the action coordinate G by the period T{G) 
or its reciprocal 1/T(G). To match the notation in Siegel-Moser, let y = 1/T(G) and x = <j). Then 
one finds that the Poincare maps take the form 



xi = x + 2vre y + e 2 f(x, y, e), 
V\ = V + e 2 g(x,y,e). 



(3.14) 



This is of the form (j3.6j) with 7 = 27re and with the perturbing functions / = e 2 f and g = e 2 g of 
order 0(e 2 ). 

We are now in a position to apply Moser's theorem to obtain: 

Theorem 3.2. Fix an annulus A of the form c < G < d with < c < d. Then for e sufficiently 
small, the Poincare map P e of the Nose-Hoover system has infinitely many invariant curves close 
to the level curves of G in A. In particular, the corresponding flow is not ergodic. 

The proof consists in checking the remaining hypotheses of Moser's theorem. Note that the 
differential equation ()3.4|) is real-analytic in (9,t, a, e). It follows that the Poincare map P e (r,a) is 
real-analytic in (r, a, e). The location of the fixed point (r e , a e ) is a real-analytic function of e (for e 
sufficiently small) and so composing the Poincare map with a translation to move this point to (1, 0) 
preserves analyticity. The action angle variables (</>, G) are independent of e and are real-analytic 
in (r, a) away from the fixed point itself. Finally, the period function T{G) is analytic in G. Thus 
the map (|3.14|) is real-analytic with respect to (x, y, e) for < y < Y, |e| < eo where Y > 0, eo > 
are constants. Since x is an angular variable, the map is also periodic with respect to x with period 
2vr. 

The annulus A contains an annulus of the form 0<a<y<6<y. If b — a = k 7^ 1, then we 
can set y = y'/k where y' is a new variable. The map will retain the same form except that 7 = 2ne 
is replaced by 7' = 2Tre/k. Once a and b are fixed, the Poincare map admits a complex-analytic 
extension to some domain of the form E = D X {|e| < ei} C C 3 where D is a closed ("(-neighborhood 
of R x [a, b] in C 2 , and 5 > 0,ei > are constants. For any fixed choice of e with |e| < ei, the 
perturbing functions / and g in (|3.14j) will have complex analytic extensions to D. We now fix the 
Diophantine constants Co and fj, in (|3.7j) and choose any e > 0. Letting 5 > be the constant in (|3.8|) 
guaranteed by Moser's theorem, the theorem shows that for any yo E [a, b] satisfying (j3.7|) . and for 
any real analytic functions /, g with complex extensions to D which satisfy (|3.8jl . the corresponding 
map admits an invariant curve e-close to y = yo- 

We now let / = e 2 f and g = e 2 g where we think of e as fixed, and we let K = sup s |/| + \g\. 
Then |/| + \g\ < e 2 K on D. Since 7 = 2vre (or 7' = 2vre/A;), it follows that (jUSJ) holds for all e 
sufficiently small. 

To complete the proof, it only remains to verify that the period satisfies T'(G) > 0. We first 
introduce a new variable a = In r. Then (|3.9jl becomes: 

a' = —a, 

(3.15) 

a' = e a - 1, V 7 

which is a planar Hamiltonian system with Hamiltonian function 

G(a, a) = \a 2 + e CT - 1 - a. 
In fact, it is of the classical kinetic plus potential form with potential 



V(a) = e a - 1 - a 
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Figure 2. Numerically computed orbits of the Poincare map of the plane 9 = 
mod 2vr for the ODE (J33J); e = 0.1 (left) and e = 1.0 (right). 

(except that in the ODE Q3.15|) . a plays the role of the momentum variable and a the role of the 
position). The equilibrium point is now at (ex, a) = (0,0). 

The behavior of the period as a function of energy for such systems is a well-studied problem. 
A result of Chicone [2] shows that T{G) will be a strictly increasing function of energy G provided 

that — t-tt^ is a strictly convex function (except at a = 0). This condition reads 
(V'(a)y 

6VV" 2 - W l2 V" - 2VV'V" > 0, 

except at a = 0. It is not hard to check that this is true for the potential V(a) above. 

Figure [2] shows numerically computed Poincare maps of the ODE (j3,3[) for two values of e. The 
Poincare map for e = 0.1 has invariant curves close to the level curves of G in Figure ^ Apparently 
the Poincare map for e = 1 still has many invariant curves although they are not particularly close 
to the level curves of G. If these really do exist, then the Nose-Hoover system is non-ergodic even 
for e = 1. 

It is interesting to remark that, for e = 1, the Poincare map invariant curves are sometimes 
composed of a set of islands (instead of being simple closed curves). This is the case for instance 
with the initial condition (9,t,q) = (0,2.42,0) which corresponds to (q,p, £) = (2.2,0,0) (7 islands 
can be seen on the right hand side of Figure EJ). Starting from this initial condition, we numerically 
integrate (|3.1j) with a time step At = 0.001 for 5.10 7 time steps with the algorithm proposed in [9]. 
This second-order algorithm is based on an operator splitting technique. It preserves the measure 
()2.3|) as well as the time-reversibility of ()3.1j) . Figure El shows the projection of the time trajectory 
of the (g,p)-plane. With this initial condition, the trajectory seems not even to sample a ring of 
the form {(q,p) £ M 2 : c < q 2 + p 2 < C} for some < c < C, but only a part of it. 

4. Harmonic oscillator coupled with a Nose-Hoover chain 

In the previous section, we have proven that a one-dimensional harmonic oscillator coupled to 
a Nose-Hoover thermostat is not an ergodic system with respect to the Gibbs measure, when the 
mass Q of the reservoir is large. As mentioned above, numerical observations of this fact have 
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q 



Figure 3. Projection on the (g,p)-plane of the numerically computed trajectory of 
1)3.1)1 with e = 1, starting from the initial condition (q,p,£) = (2.2,0,0). 



already been reported in [4,8,14], even for moderate values of Q such as Q = 1. As a consequence, 
it is known that one should be cautious when making use of the Nose-Hoover equations ([2.1)1 to 
compute phase space integrals such as 1)1.1)1 . 

One way that has been proposed to circumvent this difficulty is to generalize the Nose-Hoover 
equations to the so-called Nose- Hoover chain equations [8]. The idea consists in coupling the 
physical variables (q,p) with a first thermostat as in 1)2.1)1 and to then couple this thermostat with 
a second one, which can be coupled to a third one, and so on. The variables now include M ex t 
additional scalar variables, £j for j = 1, . . . ,M ext , where the number M ext of thermostats can be 
freely specified. The Nose-Hoover chain dynamics is given by 

dqj _ Pi 
dt mi ' 
dpi _ 
dt 

~dT^\^m i ~'^ H rQ~2^ (4-1) 



dt 



dt 




Cj for j = 2,...,M f 



cxt 



where the masses Q%, . . . , QM ext are f ree parameters that can be arbitrarily specified. 
For the Nose-Hoover system 1)4.1)1 with z = (q,p, £i, . . . ,^Af cxt )> it is easy to check that 



p(q, P, ,t 



exp —(3 



Mext c2 
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satisfies the condition (|2.2|) . so an invariant measure for the Nose-Hoover system (|4.1j) is given by 



drtsrHc(g,P,£i,--- ,6vf cx t) 



We can now calculate that 



exp ^— (3 




JLl 


^ dqdpd^i . . 


■ ^M ext 


J exp 


Mext 


2Qj 


\ dqdpd^i . 


• • ^C-Mext 



(4.2) 



r M 2 

J f-f mi 

i=i 



nMp' 1 , 



and 



/ 



< 2 

Qj 



3 dfJ-NHC = (3 1 for j = 1, 



ext • 



Hence, we can see that the evolution of £i in (|4.1|) is controlled by the difference between the 



instantaneous value of twice the kinetic energy YldLi an< ^ avera g e value with respect to the 



in, 



invariant measure cfy/NHC and the evolution of £j for j = 2, . . . , M ext in (|4.1|) is controlled by the 
difference between the instantaneous value of twice the "kinetic energy" (,j_i/Qj-i and its average 
value with respect to the invariant measure eZ/iNHC- 

If (|4.1|) is ergodic with respect to dfiNRC, computations similar to the ones performed in Section 
12 show that 

lim 1/ A(q(t),p{t))dt= [ A(q,p)dfx(q,p), (4.3) 

where (q(t), p(t), . . . ,(,M cxt (t)) is a solution of the Nose-Hoover chain dynamics and dfi(q,p) 
is the Gibbs measure (|1.2|) . Hence, the Nose-Hoover chain dynamics can be used to compute phase 
space integrals in the canonical ensemble if Q4.1JI is ergodic with respect to c^nhc- 

In this section, we numerically study the case of a one-dimensional harmonic oscillator coupled 
to a Nose-Hoover chain of length M ext = 2. We consider the case Q\ = Qi- In view of l|4.1|l . the 
dynamics reads 



P 



?2 = e K i 



-q-e pCi, 
j 2 -1- e 2 ^ 2 , 



(4.4) 



1, 



where e = = 1/VQ2- We again introduce the action-angle variables (r, 0) for the oscillator 

(see (|3.2|) ) and rescale via ay = e£j to get 



= 1 — e«i sin 8 cos 6*, 

f = — 2er«i sin 2 6*, 
di = e(2r sin 2 — 1 — 010:2), 
d 2 = e(«i - 1)- 



(4.5) 
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For small e, the corresponding averaged system reads after rescaling time, 

f = — TQl, 

di = r - 1 - aio-2, (4.6) 
a2 = af — 1. 

In the case of the Nose-Hoover equation, the averaged system is ()3.9[) . and the analysis conducted 
in Section |31 is based on the knowledge of a first integral for ()3.9|) . namely the function G(r,a) 
defined by (|3,10|) , In the case of (|4.6|) . we were not able to find such a first integral. Let us follow 
a different route and study the system by numerically computing the Poincare return map of the 
plane S^hc = {(r, 01,02) ■ 02 = 0} for the ODE (j4.6|) . Two trajectories of the return map with 
different initial conditions are shown in Figure 0J Some initial conditions lead to trajectories which 
lie on invariant curves (see right-hand side of Figure 0J and the corresponding values of t that are 
sampled are bounded from below and isolated from 0: there exists r^ =0 > such that for all t we 
have r(i) > t^ =0 . For other initial conditions, the trajectory does not seem to be confined to a 
curve (see left-hand side of Figure EJ), but it still does not sample the entire plane. 




Figure 4. Numerically computed orbits of the Poincare return map of the plane 
02 = for the ODE (j4.6j) . Left: the initial condition is (r, a\, 02) = {Qq/^, 0, 0) with 
go = 0.5; a similar trajectory is obtained for go = 1-3 and go = 1-5. Right: the initial 
condition is (r, 01,02) = (go/2,0,0) with go = 1.1; a similar trajectory is obtained 
for g = 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95 and 1.0. 

We now consider the ODE ()4.5j) and study whether the behaviour we observe in the case e — > 
still persists in the case of a small but positive e. We choose Qi = Q2 = 10, that is e = 1/VT0 = 
0.316, and we numerically integrate (|4.4j) with the second-order algorithm proposed in [9] from the 
initial condition (q,p, C11C2) = (1.1,0,0,0). This condition corresponds to an initial condition for 
which the Poincare return map of (|4.6|) seems to have invariant curves (see Figure @J right-hand 
side). Figure El shows the trace of the time trajectory of 1)4.5(1 on the plane Snhc = {(#> t, ati, 02) ■ 
02 = 0}. We can see the topological similarity between curves on Figures |U (right-hand side) 
and El although they are quantitatively quite different. Figure El shows the projection of the 
same trajectory on the (g,p)-plane. We see that the values of (q,p) that are sampled still satisfy 
r(t) = (q 2 (t)+p 2 (t))/2 > T e c for some t% > and all t > 0. Thus, the lower bound on the values of r 
that we observe for the averaged dynamics persists in the case e = 1 /VW. For the initial condition 
we consider here, t^ =0 = 0.188 whereas = 0.194. Because of this lower bound, we obtain a 
contradiction with (J4.3|) . and hence this numerical experiment indicates that a one-dimensional 
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harmonic oscillator coupled to a Nose-Hoover chain of two thermostats is not always an ergodic 
system. We observe this non-ergodicity for many different initial conditions (including all those 
listed in Figure 0J), and for many different values of Q > 10. 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 

r 

Figure 5. Trace of the trajectory of (|4.5|) with e = 1/a/10 on the plane a2 = 0. 
The initial condition is (9,r, ol\,oli ) = (O.gg/2,0,0) with g = 1-1- 

In the case Q = 1, we did not find any initial condition such that the values of r(i) are iso- 
lated from 0. Results obtained with Q = 1 and the same initial condition as previously, namely 
(<Z>P>£i> £2) = (1-1) 0,0,0), are shown on Figured We compare the theoretical distributions of the 
angular variable 6 and of the amplitude variable r = \J q 2 + p 2 (as given by the Gibbs measure 
(|1.2jl ) with the empirical distributions obtained from the time trajectory. Since we work with (3 = 1, 
the theoretical distributions respectively read 

Zt a he g o(#) = ^[0M d9 > fZo( r ) = rex P (-y) Woo) dr, 

where 1[o,2tt] is th e characteristic function of [0, 2tt]. The numerical distributions are denoted /num(#) 
and fn^S(r)- They have been computed from a trajectory of length T = 2.5 10 6 , with a time step 
of 2.5 10 -3 (the trajectory is thus composed of 10 9 time steps), by partitioning the sampled interval 
into 100 bins (a cutoff of r c = 4 has been used). Note that all distributions have been normalized 
so that their integrals are equal to 1. 



3 




-3 1 1 1 1 ^ 1 1 ^ 1 

-4 -3 -2 -1 1 234 

q 

Figure 6. Projection on the (g,p)-plane of the trajectory of Q4.5j) with e = l/\/l0, 
starting from the initial condition (9,T,a\,a2 ) = (0,^/2,0,0) with q = 1.1. 
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For the angular variable, we plot on Figure [7| both distributions f^ o (0) and We 
can see some small oscillations of the empirical distribution around the theoretical uniform dis- 
tribution. Regarding the amplitude variable, it is more convenient to directly plot the difference 
l^thco( r ) ~~ fnum( r )\- Again, we can see some oscillations, which are small (let us recall that the 
maximum value of f^^(r) is = 0.6). We have checked that the error does not change 

when the time step is reduced, the total length T of the trajectory being kept fixed. 



0.1585 - 



0.158 



i 1 1 1 1 r 

numerical _ 

theoretical 





vV + P 2 



right: difference 
The numerical 



Figure 7. Left: Numerical and theoretical distributions of 
between the numerical and theoretical distributions of r 
distribution is obtained from the simulation of Q4.5|) with e = 1, starting from the 
initial condition (q,p, £1,^2) = (1.1,0,0,0). 



To study the evolution of the error with T (the time step being now kept fixed), the indicator we 
consider is the star discrepancy. Recall that the star discrepancy Djy of a sequence x = {x n }i< n <N 
with values in [0, l] 2 is defined as [6] 



D* N (x) 



sup 

ye [0,1] 2 



1 ^ 



[0,y]t 



Volume ([0, y\) 



n=l 



(4.7) 



where, for 2-dimensional vectors y and z, we write y < z when y^ < Zi for all 1 < i < 2, and note 
that [0,y] = {z G [0, l] 2 , z < y}. The fact that D* N (x) -> when N — > 00 is equivalent [6, p. 15] to 
the fact that, for any Riemann integrable function A defined on [0, l] 2 , 



1 N f 



dx. 



In addition, for functions A which have bounded variations Vhk(^) in the sense of Hardy and 
Krause [11], the following error estimate holds true: 



1 N f 
-J2A{x n )- / A{x)dx 



< Vhk(A)D* n (x). 



(4.8) 



If A £ C 2 ([0, l] 2 ), then its variation Vh_k(A) has a simple expression: 



VkK(A) 



[0,1] 2 



d 2 A 



dx\dx2 



dx + 



dx-^ l) 



dx\ + 



■3— (M2) 
ox 2 



dx? 
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In view of l|4.8|) . we can see that the convergence of D* N {x) toward implies the convergence of 
the empirical average of A towards its statistical average, and the rate of convergence of D* N {x) 
gives information about the convergence rate of the observable average. 

Here, we work with the discrepancy criterion 

N 



D* N ({e n ,r n }) 



sup 

(6»,r)e[0,27r]x[0,r c ] 



Jf E Me] (°n) l[o,r] (r n ) ~ y_ =Q y_ =Q ^ exp 



n=l 



— ] dOdf 



(4.9) 



where the sample (9 n ,r n ) is obtained by the numerical integration of the ODE with a time step of 
At = 2.5 10~ 3 . We again set the cutoff amplitude at r c = 4. Note that the integral that appears 
in (|4.9|) can be exactly computed. In practice, D* N is approximated by considering the supremum 
only over (0, r) of the form (9 k , r{) with 6 k = 2/cvr/lOO, r, = /r c /100, and 1 < k, I < 100. 

We have considered 8 different initial conditions, giving rise to 8 different samples {0 n ,r n ), and 
for each of them, we have computed the star discrepancy for several values of N. Results are 
shown on Figure El where we plot the mean of the computed discrepancies as a function of N. A 

least mean-square fit shows that the mean discrepancy decreases as q 4g3 . 




FIGURE 8. Star discrepancy (|4.9j) as a function of the sample size N, and comparison 
to the least mean-square fit D* N « 11.1/iV 83 . 

We thus see from our numerical experiments with these thermostat masses and these initial con- 
ditions that the numerical distribution converges to the theoretical distribution, and the numerical 
results are thus consistent with ergodicity. However, more extensive numerical testing would be 
needed to assess whether this would be the case for any initial condition or any smaller thermostat 
masses. 
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